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Abstract.  The  surface  mixed  layer  of  the  ocean  is  often  characterized 
by  density  compensation  between  the  horizontal  temperature  and  salinity 
gradients.  In  this  contribution  we  present  a combination  of  theoretical 
arguments  and  numerical  simulations  to  investigate  how  compensation  might 
emerge  as  a result  of  processes  at  work  within  the  mixed  layer.  The  dynamics 
of  the  mixed  layer  are  investigated  through  a simple  model.  The  model 
consists  of  a pair  of  coupled  advection-diffusion  equations  for  heat  and  salt. 
The  coupling  arises  through  a nonlinear  diffusion  operator  proportional  to  the 
buoyancy  gradient,  which  parameterizes  the  combined  effect  of  slumping  and 
mixing  of  small-scale  horizontal  buoyancy  gradients.  Numerical  solutions  of 
the  mixed  layer  model  show  that  the  nonlinear  diffusion  creates  compensation 
between  the  temperature  and  salinity  gradients,  while  the  stirring  field 
maintains  alignment  between  the  two  gradients.  The  results  of  this  work 
suggest  a new  parameterization  of  the  horizontal  fluxes  of  heat  and  salt  for 


numerical  models  of  the  mixed  layer. 

1.  Introduction 

Observations  show  that  the  thermohaline  structure 
of  the  surface  mixed  layer  (ML)  of  the  ocean  is  largely 
compensated.  In  other  words,  temperature  and  salin- 
ity fronts  coincide  so  that  the  resulting  density  contrasts 
are  small  relative  to  the  individual  contributions  of  heat 
and  salt.  This  phenomenon  has  been  known  for  some 
time  for  certain  fronts  at  scales  of  a few  tens  to  one 
hundred  kilometers  {Roden,  1975;  Rudnick  and  Luyten, 
1996).  Recent  high-resolution  observations  have  shown 
that  compensation  exists  down  to  horizontal  scales  of 
tens  of  meters  in  the  North  Pacific  (Rudnick  and  Fer- 
rari, 1999;  Ferrari  and  Rudnick , 2000)  and  throughout 
the  global  ocean  on  scales  of  kilometers  ( Rudnick  and 
Martin,  2001).  An  example  from  a horizontal  tow  in  the 
ML  of  the  Subtropical  North  Pacific  is  given  in  Figure  1. 
Notice  how  almost  all  fluctuations  of  temperature  are 
mirrored  in  salinity  so  that  density  gradients  are  mini- 
mized. 

One  explanation  of  these  observations  is  that  atmo- 
spheric forcing  conspires  to  create  and  juxtapose  water 
masses  with  compensating  properties.  However  the  ra- 
tio of  heat  to  freshwater  density  fluxes  is  variable  in 
large  scale  maps  ( Schmitt  et  al,  1989)  and  in  time  se- 
ries at  a point  ( Weller  et  al.,  1985),  so  internal  ocean 


dynamics  is  required  to  account  for  the  observed  com- 
pensation. Young  (1994)  and  Ferrari  and  Young  (1997) 
propose  a more  satisfactory  explanation  that  relies  on 
regulating  mechanisms  at  work  in  the  ML.  These  theo- 
retical arguments  suggest  that  compensation  is  the  re- 
sult of  the  preferential  diffusion  of  horizontal  density 
gradients  which  occurs  because  of  the  combined  action 
of  unbalanced  motions  and  vertical  mixing. 

The  physical  explanation  of  the  theory  of  Young  and 
collaborators  is  as  follows.  Horizontal  gradients  of  tem- 
perature and  salinity  can  arise  in  the  ML  in  response 
to  non  homogeneous  atmospheric  forcing  and  entrain- 
ment of  thermocline  waters.  At  some  locations  temper- 
ature and  salinity  will  compensate  each  other  exactly, 
whereas  in  other  locations  temperature  and  salinity  will 
create  strong  horizontal  density  gradients.  Much  of  the 
ML  will  lie  between  these  two  extremes.  The  strong 
density  gradients  slump  under  the  action  of  gravity  and 
tend  to  restratify  the  ML.  Vertical  mixing  eventually  ar- 
rests this  unbalanced  motion  by  remixing  the  ML.  This 
mechanism  is  essentially  thermohaline  shear  dispersion, 
where  the  shear  is  driven  by  the  horizontal  density  gra- 
dient, and  the  vertical  mixing  results  from  the  variety 
of  processes  that  mix  the  ML.  On  the  other  hand,  com- 
pensated fronts  are  balanced  and  therefore  do  not  expe- 
rience shear  dispersion.  The  net  result  is  that  density 
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fronts  are  diffused,  while  compensated  fronts  persist. 

The  preferential  diffusion  of  horizontal  density  gradi- 
ents can  be  represented  with  mixing  parameterizations 
in  which  the  transport  of  heat  and  salt  depends  nonlin- 
early  on  the  density  gradient,  e.g.,  with  diffusivities  pro- 
portional to  some  power  of  the  density  gradient  ( Young , 
1994;  Ferrari  and  Young , 1997).  In  this  paper  we  ex- 
amine the  establishment  of  thermohaline  compensation 
by  implementing  these  nonlinear  diffusive  parameteri- 
zations in  a simple  model  of  the  ML.  Numerical  solu- 
tions show  that  the  thermohaline  structure  of  the  ML 
is  generated  by  a balance  between  the  mesoscale  strain- 
ing field,  that  acts  to  increase  temperature  and  salinity 
gradients,  and  the  nonlinear  diffusion,  that  arrests  the 
formation  of  density  gradients  but  not  of  compensated 
gradients.  In  agreement  with  observations,  temperature 
and  salinity  gradients  tend  to  be  aligned,  because  both 
heat  and  salt  are  advected  by  the  same  straining  field, 
and  compensated. 

The  mechanism  of  compensation  described  above  im- 
plicate vertically  sheared  currents  within  the  ML  and  it 
is  not  included  in  numerical  models  with  bulk  MLs.  In 
the  last  section  of  the  paper  we  show  how  to  simplify 
the  nonlinear  diffusive  parameterization  so  that  it  can 
be  implemented  in  ocean  circulation  models  to  improve 
the  representation  of  ML  thermohaline  dynamics. 

The  paper  is  organized  as  follows.  In  section  2,  we 
revisit  the  arguments  of  Young  and  collaborators  in  the 
context  of  the  parameterization  of  diapycnal  fluxes  in 
the  ML.  In  section  3,  we  describe  numerical  simulations 
used  to  test  the  nonlinear  diffusive  parameterization  of 
heat  and  salt  transports  in  the  ML.  In  section  4,  we 
suggest  a simplified  version  of  the  nonlinear  diffusive 
parameterization  to  be  implemented  in  bulk  ML  mod- 
els. Finally,  conclusions  are  offered  in  section  5. 

2.  Horizontal  transport  of  heat  and  salt 
in  the  mixed  layer 

Let  us  consider  the  dispersion  of  some  tracer  of  con- 
centration 6 in  the  ML.  We  model  the  ML  as  a vig- 
orously mixed,  shallow  layer,  characterized  by  a small 
aspect  ratio,  i.e.,  with  a depth  H much  less  than  the 
horizontal  scale.  The  main  point  here  is  that  there  are 
two  very  different  time  scales:  a fast  time  scale  ry  over 
which  the  layer  is  mixed  vertically  over  the  depth  H 
and  a longer  time  scale  t#  associated  with  horizontal 
transports. 

The  mathematical  model  for  the  transport  of  a tracer 
6 stirred  by  an  incompressible  velocity  field  u is  the 
familiar  advection-diffusion  equation, 

(1) 


together  with  appropriate  boundary  conditions.  The 
operator  on  the  RHS  represents  the  diffusion  of  tracer 
fluctuations  by  molecular  motions  and  k is  the  molec- 
ular diffusivity.  The  equation  in  (1)  is  appropriate  to 
describe  the  transport  of  6 at  scales  from  a few  millime- 
ters to  thousands  of  kilometers.  However,  the  resulting 
description  is  overly  complicated.  Our  goal  is  to  derive 
a simpler  model  that  describes  transports  in  the  ML  at 
large  scales  and  long  times  by  averaging  the  equation 
in  (1)  over  short  times  and  short  scales.  The  key  step 
in  the  analysis  is  to  find  appropriate  scales  for  the  av- 
eraging so  that  we  can  derive  a closed  equation  for  the 
averaged  concentration  6 by  folding  all  the  details  of 
the  small  scale  motions  in  a suitable  operator  V that 
depends  on  averaged  variables,  i.e., 

5*0  + u-V0  = X>(u,0),  (2) 

where  u is  the  averaged  velocity  field. 

Some  very  popular  ML  models,  referred  to  as  bulk 
models  (e.g.  Kraus  and  Turner,  1967),  choose  to  aver- 
age over  the  depth  H of  the  ML  and  the  characteristic 
time  of  vertical  mixing  ry . This  choice  is  quite  natural, 
because  the  turbulent  fluxes  that  homogenize  vertically 
the  ML  are  due  to  processes  such  as  convection  and 
Langmuir  cells,  characterized  by  coherent  eddies  which 
span  the  depth  H and  have  an  aspect  ratio  close  to  one. 
In  these  models,  the  operator  V parameterizes  all  the 
processes  that  maintain  the  ML  well  mixed  in  the  ver- 
tical. A problem  arises  when  bulk  ML  are  implemented 
in  circulation  models  that  resolve  horizontal  scales  that 
are  orders  of  magnitude  larger  than  H.  In  this  case 
one  has  to  average  the  tracer  equation  over  H in  the 
vertical,  but  over  a scale  L > H in  the  horizontal.  A 
typical  solution  is  to  parameterize  in  series  the  motions 
on  scales  shorter  than  H and  those  on  scales  between 
H and  Ll.  That  is,  the  same  operator  V is  retained 
to  describe  the  fluxes  that  mix  vertically  the  ML,  but 
a lateral  effective  eddy  diffusivity  is  introduced  to  pa- 
rameterize the  fluxes  at  larger  scales.  Here  we  show 
that  imbalanced  horizontal  motions  with  characteristic 
scales  between  H and  L act  in  parallel  with  the  turbu- 
lent motions  on  scales  shorter  than  H.  Therefore  it  is 
necessary  to  modify  the  operator  T>  and  parameterize 
all  unresolved  motions  together. 

Let  us  average  equation  (1)  over  the  depth  H in  the 
vertical,  over  a scale  L > H in  the  horizontal  and  over 
a time  th  > ry.  The  scales  L and  th  have  only  lower 
bounds,  but  are  not  specified  for  the  moment.  We  ob- 
tain the  Reynolds’  averaged  equation, 

dt&  + u-  VH0  = -VH  ■urF  + K'V%d  + T.  (3) 


1See  Garrett  (2001)  for  a discussion  of  parameterizations  of 
unresolved  motions  in  parallel  and  in  series. 


dto  + u • ve  = kv26, 


THE  TEMPERATURE-SALINITY  RELATIONSHIP  OF  THE  MIXED  LAYER 


97 


28 


28.5 


29 

Latitude  (degN) 


29.5 


36 


35.5 


35 


D 

CO 

CL 

CO 


34.5 


30 


Figure  1.  Potential  temperature  (red  line),  salinity  (blue  line)  and  potential  density  (black  line)  from  a horizontal  SeaSoar 
tow  at  50  m in  Subtropical  North  Pacific,  at  140  degrees  west,  between  28  and  30  degrees  north.  This  depth  is  in  the  middle 
of  the  local  mixed  layer.  The  vertical  axis  for  temperature  and  salinity  are  scaled  by  the  respective  expansion  coefficients  so 
that  excursions  of  temperature  and  salinity  show  the  change  they  imply  on  density. 


Here  6 and  u are  the  averaged  velocity  and  averaged 
tracer  concentration  and  6'  and  u'  are  departures  from 
those  averages.  T represents  the  flux  of  tracer  induced 
by  the  boundary  conditions  at  the  top  and  bottom  of 
the  ML.  The  notation  V#  is  used  to  remind  that  deriva- 
tives are  taken  only  in  the  horizontal,  because  the  av- 
eraged quantities  do  not  depend  on  the  vertical  coor- 
dinate. In  order  to  simplify  the  discussion,  we  assume 
that  H is  a constant  independent  of  position  (for  more 
on  this  point  see  Young,  1994;  Garrett  and  Tandon, 
1997). 

The  next  step  is  to  express  the  first  term  on  the 
RHS  (called  “eddy  flux  divergence”)  in  terms  of  aver- 
aged quantities.  Mixing-length  theories  are  a common 
way  to  achieve  this  goal.  The  argument  goes  that  a 
fluid  particle  carries  the  value  of  a conserved,  and  hence 
transferable,  tracer  for  some  length  I',  before  it  is  mixed 
with  its  new  surroundings.  We  give  a vectorial  nature 
to  1'  to  allow  for  situations  which  are  not  isotropic.  If 
the  particle  has  a concentration  of  scalar  typical  of  its 
surroundings  then  the  eddy  flux  of  tracer  6 is  given  by 

u^-uT-  XH0,  (4) 

where  it  is  assumed  that  Vh®  varies  little  over  distances 
comparable  with  the  mixing  length  1'.  The  tensor  u'l' 
defines  the  eddy  diffusivity. 

In  the  special  case  when  the  statistics  of  the  velocity 
field  are  homogeneous  and  isotropic,  the  eddy  diffusivity 
tensor  is  a constant,  and  the  eddy  transport  assumes  the 


form  of  a down-gradient  Fickian  diffusion: 

= -kXH6,  (5) 

This  kind  of  closure  is  commonly  applied  to  ML  models 
and  the  two  relevant  scalars  (temperature  and  salinity) 
are  diffused  with  the  same  eddy  diffusion  coefficient, 
and  are  independent  from  each  other. 

However  in  the  ML  there  are  lateral  inhomogeneities 
in  the  buoyancy2  field  at  scales  larger  than  H.  Hor- 
izontal buoyancy  gradients  slump  under  the  action  of 
gravity  and  drive  horizontal  eddy  fluxes.  Therefore  we 
expect  the  transport  of  tracer  to  be  in  the  direction  of 
and  to  increase  with  Vj ?J3.  This  breaks  the  assump- 
tions of  homogeneity  and  isotropy.  Therefore  a down- 
gradient  Fickian  diffusion  cannot  be  used  to  model  the 
ML  at  scales  larger  then  H.  A more  appropriate  ex- 
pression for  the  diffusivity  tensor  is, 

= 7/(|Vtf S|)VhB  vhB,  (6) 

where  7 is  a constant  and  /(|Vij.B|)  a non-dimensional 
function  whose  form  depends  on  the  details  of  the 
hydrodynamic  instabilities  that  dominate  in  the  eddy 
field.  The  expression  (6)  is  rationalized  as  follows.  Ac- 
cording to  mixing-length  theories,  the  diffusivity  tensor 

2 Buoyancy  B is  defined  as  p = po  \l  - , where  p is 

the  fluid’s  density,  po  is  a constant  reference  density  and  g is  the 
acceleration  of  gravity.  With  this  definition,  B has  the  dimensions 
of  acceleration. 
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can  be  expressed  in  terms  of  the  characteristic  velocity 
veddy  and  length  leddy  of  the  transfer  process,  that  is 
u'l'  oc  ueddy  1 eddy  In  our  case  the  length  scale  is  given 
by  leddy  = ^eddyTy,  where  Ty  is  the  time  for  which 
the  slumping  process  acts,  before  it  is  arrested  by  the 
turbulent  fluxes  that  mix  vertically  the  ML.  The  eddy 
velocity  is  in  the  direction  of  Vi/B  with  a magnitude 
proportional  to  |VtfB|.  Thus  we  get  the  expression  in 
(6),  where  the  tensor  V # B V#  B arises  from  the  direc- 
tion of  the  eddy  velocity  field  and  7/(|  V# B|)  is  a pos- 
itive semidefinite  term  that  determines  the  magnitude 
of  the  flux.  Notice  that  both  the  unbalanced  motions  at 
scales  larger  than  H and  the  turbulent  fluxes  at  scales 
shorter  than  H enter  in  the  closure  in  6.  That  is  the 
processes  of  slumping  and  mixing  act  in  parallel. 

Plugging  (6)  into  (4)  gives  the  eddy  tracer  flux, 

= -7/(|V*B|)  (VhB  • Vh6)  V* B.  (7) 

Notice  that,  even  though  the  flux  is  in  the  direction  of 
VB,  u'9'  ■ V#  < 0.  Thus  the  flux  of  tracer  tends  to 
be  down  the  tracer  gradient,  but  only  the  projection  of 
the  tracer  gradient  along  the  direction  of  the  buoyancy 
gradient  contributes  to  the  flux. 

We  now  apply  the  closure  in  (7)  to  the  advection- 
diffusion  equations  for  heat  and  salt  in  the  ML, 

dtT  + u ■ VT  = 

= 7V-[/(|VB|)(VB-Vr)VB]+JFT,  (8) 

dtS  + 1 1 • V5  = 

= 7V-[/(|VB|)(VB-VB)VB]+JTs,  (9) 

where  Ft  and  Fs  represent  the  thermohaline  fluxes 
from  the  top  and  bottom  of  the  ML.  We  dropped  over- 
bars  and  we  replaced  V#  with  V for  convenience,  but 
keep  in  mind  that  all  variables  are  averaged  over  scales 
larger  than  H and  times  longer  than  ry  and  that  deriva- 
tives are  taken  only  in  the  horizontal.  We  assume  a lin- 
ear equation  of  state  and  measure  T and  S in  buoyancy 
units  units,  so  that, 

B = T - S.  (10) 

The  nonlinear  advection-diffusion  equations  (8)  and  (9), 
together  with  (10),  form  a closed  system  whose  solution 
is  fully  determined  once  the  forcings  Ft  and  Fs  and  the 
large  scale  velocity  field  u are  prescribed. 

By  adding  and  subtracting  (8)  and  (9),  we  obtain 
closed  equations  for  buoyancy  and  spice  V = T + S 
( Veronis , 1972;  Munk,  1981),  viz., 

dtB  + u • VB  = 

= 7V  • [/(|VB|)  (VB  • VB)  VB]  + ^s,(ll) 
dtV  + u • W = 

= 7V  • [/(|VB|)  (W  • VB) VB]  + Fv,(  12) 


where  Fb  = Ft  — Fs  and  Fy  — Ft  + Fs-  We  can 
now  see  how  the  equations  in  (11)  and  (12)  model  the 
development  of  compensation  in  the  ML.  The  prod- 
uct /(|VB|)(VB  • VB)  must  be  an  increasing  func- 
tion of  | VB  I to  be  consistent  with  our  assumption  that 
eddy  fluxes  are  driven  by  buoyancy  gradients.  Under 
this  constraint,  the  nonlinear  diffusion  always  dissipates 
buoyancy,  even  more  so  when  |VB|  is  large.  Also  spice 
is  dissipated  where  |VB|  is  large.  However  large  val- 
ues of  |VVj  can  survive  in  regions  where  |VB|  is  small. 
In  terms  of  temperature  and  salinity  this  means  that 
compensated  fronts,  for  which  VT  « VB  persist,  while 
buoyancy  fronts  are  short  lived. 

Young  (1994)  and  Ferrari  and  Young  (1997)  derive 
formally  equations  of  the  form  of  those  in  (8)  and  (9) 
to  parameterize  the  transport  of  heat  and  salt  on  hori- 
zontal scales  of  a few  kilometers  in  the  ML.  These  theo- 
retical works  are  examples  of  the  closures  we  have  been 
discussing  when  the  averaging  is  done  over  the  depth  of 
the  ML,  over  horizontal  scales  of  a few  kilometers  and 
time  scales  of  a few  hours.  Nonlinear  diffusion  arises 
because  the  horizontal  transport  of  heat  and  salt  is  by 
shear  dispersion,  and  the  shear  flow  doing  the  disper- 
sion is  driven  by  slumping  horizontal  buoyancy  gradi- 
ents. The  strength  of  the  shear  dispersion  increases 
as  the  horizontal  buoyancy  gradient  squared,  that  is 
/(|VB|)  = 1 in  (8)  and  (9). 

At  scales  larger  than  the  Rossby  radius  of  deforma- 
tion Ro , unbalanced  motions  are  influenced  by  rotation 
in  the  form  of  baroclinic  instability.  Therefore,  if  one 
is  to  parameterize  the  transport  of  heat  and  salt  on 
horizontal  scales  larger  than  Ro,  say  10  km  for  a typi- 
cal ML,  the  closure  must  include  the  transports  due  to 
eddies  generated  at  baroclinically  unstable  gradients. 
Green  (1970)  and  Stone  (1972)  derived  expressions  for 
the  tracer  fluxes  generated  by  baroclinic  waves.  Their 
results  predict  that  the  baroclinic  eddy  fluxes  across 
a buoyancy  gradient  are  proportional  to  the  absolute 
value  of  the  diapycnal  buoyancy  gradient.  Green  and 
Stone  considered  only  zonally-averaged  models  and  did 
not  investigate  the  direction  of  the  fluxes.  If  their  ar- 
guments are  extended  to  two  horizontal  dimensions  to 
parameterize  diapycnal  fluxes  of  heat  and  salt  in  the 
ML,  one  obtains  nonlinear  diffusion  equations  of  the 
form  in  (8)  and  (9)  with  /(|VB|)  = |VB|_1.  Notice, 
however,  that  a full  parameterization  of  baroclinic  in- 
stability should  include  the  eddy  fluxes  along  isopycnals 
as  well  ( Marshall  and  Shutts,  1981).  This  issue  is  not 
pursued  further  here,  because  we  focus  on  the  role  of  di- 
apycnal fluxes  on  the  establishment  of  the  temperature- 
salinity  relationship  in  the  ML. 

Chris  Garrett,  during  the  meeting,  suggested  that 
symmetric  instability  might  also  drive  thermohaline 
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fluxes  in  the  ML.  Haine  and  Marshall  (1997)  used  nu- 
merical simulations  to  study  what  hydrodynamical  in- 
stabilities control  the  transfer  of  buoyancy  through  the 
ML  on  scales  of  some  tens  of  kilometers.  Their  conclu- 
sion is  that  nonhydrostatic  baroclinic  instability  pro- 
vides the  dominant  mode  of  lateral  buoyancy  transfer. 
However  symmetric  instability  plays  an  important  role 
during  the  slumping  process  by  setting  to  zero  potential 
vorticity  along  isopycnal  surfaces.  Clearly  more  work 
need  to  be  done  to  formulate  a closure  that  takes  into 
account  the  effects  of  both  symmetric  and  baroclinic 
instabilities. 

3.  Thermohaline  alignment  and 
compensation  in  the  mixed  layer 

The  nonlinear  advection-diffusion  equations  in  (8) 
and  (9)  are  now  used  to  investigate  how  compensation 
appears  in  the  ML.  Suppose  that  spatial  variations  in 
temperature  and  salinity  are  created  by  surface  fluxes 
that  vary  on  large  horizontal  scales.  Mesoscale  stirring 
will  create  small-scale  temperature  and  salinity  gradi- 
ents by  stretching  and  folding  the  large  scale  thermo- 
haline patterns.  Large  density  gradients  will  disappear 
quickly  as  a result  of  nonlinear  diffusion,  while  com- 
pensated gradients  will  persist  for  longer  times.  Thus 
we  expect  that  the  temperature  and  salinity  gradients 
present  at  small  scales  at  any  particular  moment  will  be 
typically  compensated.  We  are  now  going  to  test  this 
scenario  with  a numerical  model. 

3.1.  Numerical  model 

The  parameterization  in  (8)  and  (9)  is  tested  with  nu- 
merical simulations  in  which  temperature  and  salinity 
are  advected  using  a velocity  field  generated  by  solving 
the  equivalent  barotropic  equations  in  the  streamfunction- 
vorticity  formulation, . 

dtC  + j(4>,  0 = + "V6C  + (13) 

where  ip  is  the  streamfunction,  ( — V2ip  the  relative 
vorticity  and  J the  Jacobian  operator.  The  forcing  is 
applied  in  spectral  space  at  a scale  of  6 km  with  constant 
amplitude  and  random  phases.  The  bottom  drag  coef- 
ficient is  set  to  /i  = 3 • 10~6  s-1  and  the  hyper-viscosity 
to  v = 3 • 106  m6  s-1.  The  result  is  a two-dimensional 
turbulent  field  with  meandering  vortices  of  a diameter 
of  approximately  3 km  (half  the  forcing  scale)  and  RMS 
velocities  of  0.1  m s-1  (Figure  2).  The  domain  of  inte- 
gration is  a biperiodic  square  of  51.2  x 51.2  km2  with  a 
mesh  of  100  m.  This  is  a poor  model  of  the  mesoscale 
dynamics  of  the  ML.  In  particular  we  are  neglecting 
feedbacks  between  the  buoyancy  and  the  velocity  fields. 
But  our  goal  is  to  show  that  compensation  develops  at 


Figure  2.  Snapshot  of  the  vorticity  field  obtained  by  in- 
tegrating the  equivalent  barotropic  equations.  The  typical 
size  of  vortices  is  about  3 km,  that  is  half  the  wavelength  of 
6 km  at  which  the  vorticity  equation  is  forced. 


scales  small  scales,  regardless  of  the  details  of  the  stir- 
ring field  and,  in  this  context,  the  model  in  13  suffices. 

The  temperature  and  salinity  equations  in  (8)  and 
(9)  are  integrated  with  /(|VB|)  = 1,  that  is  we  use  the 
closure  in  Young  (1994)  and  Ferrari  and  Young  (1997). 
The  value  of  7 is  set  to  10 14  m2s3  appropriate  for  typical 
ML  parameters  (details  in  Ferrari  and  Young,  1997). 
However  the  qualitative  results  discussed  in  the  rest  of 
this  paper  do  not  depend  on  the  particular  choice  of 
/(|V2?|). 

Temperature  and  salinity  are  forced  with  orthogonal 
sinusoidal  patterns,  that  is  we  set  Ft  = Fq  cos  qx  in  the 
RHS  of  (8)  and  Fs  = Fq  sin  qy  in  the  RHS  of  (9).  The 
sinusoids  have  a wavelength  equal  to  the  domain  size, 
i.e.,  q = 27r/51.2  km-1.  The  amplitude  Fo  is  chosen 
such  as  to  have  thermohaline  fluctuations  of  1°C  and 
0.35  psu,  at  the  scale  of  the  domain.  These  forcings  do 
not  to  impose  any  correlation  between  temperature  and 
salinity  fluctuations.  Further  details  on  the  numerical 
code  are  given  in  Ferrari  and  Paparella  (2001). 


3.2.  Complex  density  ratio 


It  is  common  to  quantify  compensation  in  terms  of 
the  density  ratio,  defined  as  the  change  in  buoyancy  due 
to  temperature  divided  by  the  change  in  buoyancy  due 
to  salinity, 


Rid  = 


bVT 

i • vs’ 


(14) 


where  temperature  and  salinity  are  defined  in  buoyancy 
units,  and  l is  the  direction  along  which  the  cut  is  taken. 
In  two-dimensions,  it  is  convenient  to  introduce  a 
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complex  density  ratio  as, 

d _ Tx  + iTy 

SX+iSy 


(15) 


The  complex  density  ratio  has  both  a magnitude  and  a 
phase,  R = |1?|  exp(i(?i):  |1?|  is  the  ratio  of  the  magni- 
tudes of  the  temperature  and  salinity  gradients  and  <j> 
is  the  angle  between  them.  If  the  gradients  are  parallel 
(<f>  = 0°)  or  antiparallel  (<p  — 180°),  there  is  thermoha- 
line alignment  and  the  definition  of  R in  the  complex 
plane  is  equivalent  to  that  in  (14)  regardless  of  the  ori- 
entation of  l.  When  |i?|  > 1,  the  change  in  buoyancy 
due  to  temperature  is  greater  than  the  change  in  buoy- 
ancy due  to  salinity  along  the  direction  aiVB,  that 
is  | VB  ■ VT|  > | VB  ■ V5|  and  the  buoyancy-front  is 
temperature-dominated.  The  opposite  is  true  if  |i?|  < 1 
and  the  buoyancy-front  is  salinity-dominated.  The  par- 
ticular case  |1?|  = 1 and  4>  — 0°  describes  thermohaline 
compensation. 

Because  the  magnitude  of  the  complex  density  ratio 
is  infinite  when  the  salinity  gradient  vanishes,  we  char- 
acterize fronts  in  terms  of  the  phase  <f>  and  the  Turner 
angle, 

Tu  = arctan  |i£|,  (16) 

choosing  the  branch  where  0 < Tu  < 7r/2.  All  statistics 
will  be  computed  in  terms  of  (j>  and  Tu.  For  convenience, 
results  are  discussed  in  terms  of  <t>  and  | J?|,  because  their 
values  are  more  familiar. 

In  the  following  we  use  the  joint  pdf  P(Tu,  <f>)  to  de- 
scribe the  degree  of  alignment  and  compensation  in  the 
ML.  The  joint  pdf  is  normalized  according  to 

rff/2  r2n 

/ dTu  / &<p  P(Tu, <{>)  = 1.  (17) 

Jo  Jo 


3.3.  Results  of  numerical  simulations 

For  the  simulations  we  use  kilometers  to  measure  dis- 
tances and  hours  to  measure  time.  Therefore  vorticity 
is  given  in  h_1  and  buoyancy  in  km  h-2.  We  set  to  zero 
the  initial  vorticity,  temperature  and  salinity.  After  an 
initial  transient  of  several  eddy  turnover  times,  kinetic 
energy,  enstrophy,  temperature  and  salinity  variances 
settle  to  a constant  value;  i.e.  the  system  reaches  an 
equilibrium  between  the  variance  input  by  the  forcing 
at  large  scales  and  dissipation  at  small  scales. 

In  Figures  3 and  4,  we  show  snapshots  of  spice  and 
buoyancy  700  h after  the  beginning  of  the  simulation. 
It  is  difficult  to  recognize  in  these  snapshots  the  large 
scale  patterns  of  buoyancy  and  spice  imposed  by  the 
forcing  described  in  section  3.1.  But  the  sinusoidal  pat- 
terns emerge  clearly  if  one  averages  the  fields  over  times 
of  the  order  of  a few  hundred  hours.  At  small  scales  the 


two  fields  are  remarkably  different.  A comparison  of  the 
black  contours  in  the  two  figures  shows  that  gradients 
of  spice  are  sharper  than  those  of  buoyancy:  buoyancy 
contours  are  evenly  spaced,  while  spice  contours  are  ex- 
tremely packed  in  a few  regions  and  widely  spaced  in 
others.  Sharp  spice  gradients  with  no  signature  in  buoy- 
ancy imply  VT  « VS,  i.e.,  thermohaline  compensation. 

The  small  scale  variability  in  Figures  3 and  4 is 
produced  by  stirring  the  large  scale  thermohaline  pat- 
terns. The  temperature  and  salinity  gradients,  com- 
puted across  the  grid  spacing  of  100  m,  are  typically 
aligned.  Alignment  occurs  because  the  isolines  of  T 
and  S are  stretched  by  the  same  stirring  field  and  thus 
both  tracers  end  up  with  gradients  pointing  in  the  same 
directions  ( Hua , 2001).  This  is  shown  through  the  joint 
pdf  'P(Tu,  tp)  (Figure  5).  The  overwhelming  majority  of 
points  in  the  pdf  have  angles  very  close  to  either  (p  — 0° 
or  <f>  = 180°.  But  this  is  not  the  whole  story,  because 
not  all  values  of  |i?|  are  equally  probable  along  those 
two  angles.  The  pdf  has  a clear  peak  at  R = 1.  This 
is  the  signature  of  nonlinear  diffusion  which  selectively 
dissipates  all  gradients  whose  density  ratio  is  different 
from  one  and  establishes  compensation.  Stirring  alone 
does  not  produce  a single  peak  in  the  pdf,  because  it 
acts  only  on  the  relative  orientation  of  VT  and  VS  but 
not  on  the  ratio  of  their  magnitudes.  This  was  checked 
by  running  a simulation  in  which  the  nonlinear  diffusion 
was  set  to  zero.  In  this  limit,  the  pdf  P(Tu,  </>)  is  indeed 
collapsed  along  the  angles  <p  — 0°  and  (j>  = 180°,  but  it 
does  not  have  a single  mode. 

Compensation  is  not  maintained  always  and  every- 
where in  the  domain.  There  are  regions,  in  Figures  3 
and  4,  where  buoyancy  and  spice  gradients  are  compa- 
rable. This  happens  when  the  stirring  field  momentarily 
creates  large  buoyancy  gradients  at  small  scales.  These 
gradients  do  not  persist  for  long,  though,  because  non- 
linear diffusion  restores  compensation  in  a few  hours. 
At  any  time,  a one  dimensional  cut  through  the  do- 
main shows  many  compensated  fronts  and  some  rare 
buoyancy  fronts.  This  result  agrees  with  the  thermo- 
haline structure  found  by  Ferrari  and  Rudnick  (2000) 
in  the  ML  of  the  Subtropical  North  Pacific  (Figure  1), 
where  almost  all  temperature  and  salinity  fluctuations 
are  compensated. 

4.  Implications  for  numerical  models  of 
the  mixed  layer 

In  the  previous  sections  we  have  suggested  that  the 
thermohaline  compensation  observed  in  the  ML  is  con- 
sistent with  preferential  diffusion  of  horizontal  buoy- 
ancy gradients.  The  theoretical  argument  implicates 
vertically  sheared  currents  within  the  ML  as  the  agent 
which  produces  the  preferential  horizontal  transport  of 


THE  TEMPERATURE-SALINITY  RELATIONSHIP  OF  THE  MIXED  LAYER 


101 


Figure  3.  Snapshots  of  spice  at  the  same  time  of  figures  2 
and  4.  The  colored  pattern  hints  at  the  large-scale  sinu- 
soidal checkerboard,  imposed  on  the  spice  field  through  the 
thermohaline  forcing. 
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Figure  5.  Joint  pdf  ^(Tu,  <j>)  of  the  complex  density  ratio  R 
defined  in  section  3.2.  The  azimuthal  position  in  the  radar 
plot  indicates  the  angle  between  the  temperature  and  the 
salinity  gradients,  while  the  radial  displacement  is  the  ratio 
of  their  magnitudes.  The  pdf  is  an  average  over  200  hours 
during  which  the  simulation  was  in  equilibrium.  Nearly  all 
points  lie  along  the  line  of  alignment,  that  is  at  angles  <j>  = 0° 
and  <f>  = 180°.  The  maximum  of  the  pdf  is  at  R — 1 and 
shows  that  thermohaline  fronts  are  typically  compensated. 
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Figure  4.  Snapshots  of  buoyancy  at  the  same  time  of 
figures  2 and  3.  The  colored  pattern  hints  at  the  large- 
scale  sinusoidal  checkerboard,  imposed  on  the  buoyancy  field 
through  the  thermohaline  forcing.  At  small  scales  buoy- 
ancy shows  less  structure  and  milder  gradients  than  spice. 
Regions  with  large  spice  fluctuations  with  no  signature  in 
buoyancy  are  the  trademark  of  compensation. 


density.  Numerical  models  with  bulk  MLs  do  not  in- 
clude this  physics.  The  purpose  of  this  last  section  is 
to  suggest  a simple  parameterization  of  ML  horizon- 
tal transports  that  might  improve  the  representation  of 
thermohaline  dynamics  in  numerical  models  with  bulk 
ML. 

Bulk  ML  parameterizations  ignore  the  potential  en- 
ergy stored  in  horizontal  buoyancy  gradients.  But  we 
contend  that  the  release  of  this  potential  energy  plays 
an  important  role  in  establishing  compensation.  The 
system  in  (11)  and  (12)  describes  the  horizontal  dynam- 
ics of  the  vertically-averaged  fields  and  could  be  imple- 
mented in  a bulk  ML  model.  However  the  nonlinear 
diffusive  terms  on  the  RHS  of  (11)  and  (12)  are  difficult 
to  integrate  numerically.  In  small  regions  with  large 
buoyancy  gradients  the  diffusive  constraint  on  the  time 
stepping  becomes  severe  ( Ferrari  and  Paparella,  2001) 
and  the  whole  calculation  proceeds  very  slowly.  The 
path  we  follow  here  is  to  derive  a substitute  model  that 
retains  the  basic  physics  of  (11)  and  (12),  but  that  is 
easy  to  integrate  numerically.  That  is  we  write  a linear 
model  that  diffuses  horizontal  buoyancy  gradients  more 
efficiently  than  spice  gradients  in  the  following  way, 


= kbV2B, 

- kvV2V. 


Bt  + u • VB 
Vt  + u • VV 


(18) 

(19) 
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By  setting  kb  Ky  buoyancy  gradients  decay  faster 
than  spice  gradients.  The  main  advantage  of  the  model 
in  (18)  and  (19)  over  that  in  (11)  and  (12)  is  that  the 
diffusivities  are  independent  of  the  buoyancy  gradients 
and  therefore  the  diffusive  constraints  on  the  time  step- 
ping do  not  grow  unbounded. 

In  order  to  match  observations,  the  two  diffusivi- 
ties kb  and  Ky  must  be  chosen  such  that  compensa- 
tion happens  mostly  at  scales  below  10  km  ( Ferrari  and 
Rudnick,  2000).  That  is  the  dissipation  cutoff  scale  for 
buoyancy  must  be  of  the  order  of  10  km,  while  the  dis- 
sipation cutoff  scale  for  spice  must  be  smaller.  The 
dissipation  cutoff  scale  for  buoyancy  can  be  estimated 
as, 


where  a is  the  RMS  strain  rate  of  to  the  mesoscale  eddy 
field  u.  A reasonable  strain  rate  in  the  ML  is  of  the  or- 
der of  10"5  s-1.  By  imposing  Ldus  ~ 10  km,  it  follows 
that  kb  » 103  m2  s-1.  The  choice  of  kv  is  somewhat 
arbitrary,  but  it  should  be  a couple  of  orders  of  magni- 
tude smaller  than  kb,  so  that  there  is  at  least  a decade 
between  the  cutoff  scales  of  spice  and  buoyancy. 

The  final  step  is  to  write  the  linear  model  in  (18)  and 
(19)  in  terms  of  temperature  and  salinity,  by  using  once 
more  the  linear  expressions  for  buoyancy  B — T—S  and 
spice  V = T + S, 

Tt  + n-VT  = k+V2T  — k_V2S,  (21) 

St+n-VS  = k+V2S  — k_V2T,  (22) 

where  k+  = ( kb  + Ky)/2  and  k_  = (kb  — Ky)/2.  The 
coupling  between  the  salt  and  heat  fluxes  in  (21)  and 
(22)  is  formally  similar  to  the  Soret  and  DuFour  effects 
that  operate  on  a molecular  level  (Caldwell,  1973).  The 
main  difference  is  that,  in  the  present  case,  all  terms 
in  the  RHS  of  (21)  and  (22)  are  of  the  same  order  and 
none  can  be  neglected.  Only  for  kb  = Ky,  the  coupling 
between  temperature  and  salinity  disappears. 

The  parameterization  in  (21)  and  (22)  has  been 
tested  by  integrating  the  equations  with  the  velocity 
field  and  thermohaline  forcings  described  in  section  3.1. 
The  lower  panel  of  Figure  6 shows  the  pdf  of  the  merid- 
ional density  ratio  at  a scale  of  3 km  obtained  with  the 
linear  model.  The  pdf  has  a clear  peak  at  R = 1,  as 
in  the  observations  (upper  panel  of  Figure  6).  Notice 
that  the  large  scale  density  ratio  for  the  same  simula- 
tion is  uniform.  Thus  compensation  is  a result  of  the 
parameterization  in  (21)  and  (22),  and  is  not  due  to  the 
external  forcing. 
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Figure  6.  Probability  density  function  of  the  horizontal 
mixed-layer  Turner  angle  across  a distance  of  3 km  along 
25-35  degrees  in  the  North  Pacific  (upper  panel)  and  across 
the  same  distance  in  a simulation  with  the  proposed  mixed 
layer  parameterization  (lower  panel).  The  values  of  Turner 
angle  and  density  ratio  are  indicated  on  the  upper  and  lower 
axes.  The  pdfs  have  a peak  close  to  R = 1 and  represent 
thermohaline  fields  characterized  by  buoyancy  compensa- 
tion. 

5.  Conclusions 

We  have  shown  that  the  ubiquitous  compensation 
of  thermohaline  gradients  observed  in  the  ML  is  consis- 
tent with  the  theoretical  arguments  of  Young  (1994)  and 
Ferrari  and  Young  (1997).  Compensation  can  be  ex- 
plained as  the  preferential  diffusion  of  horizontal  buoy- 
ancy gradients  which  occurs  because  unbalanced  motion 
due  to  these  gradients  is  stronger  in  the  ML  than  in  the 
more  nearly  geostrophic  interior.  The  horizontal  pres- 
sure gradients  associated  with  the  buoyancy  gradients 
produce  “exchange  flows”  which  act  to  restratify  the 
ML  in  the  vertical.  The  turbulent  fluxes,  that  contin- 
uously mix  the  ML  in  the  vertical,  oppose  the  restrati- 
fication and  weaken  the  horizontal  buoyancy  gradients. 
This  process  is  essentially  shear  dispersion  of  buoyancy, 
where  the  shear  flow  is  driven  by  the  density  gradients 
themselves. 

The  theoretical  arguments  of  Young  and  collabora- 
tors implicate  that  eddy  fluxes  of  heat  and  salt  in  the 
ML  are  in  the  direction  of  the  buoyancy  gradients  and 
act  to  weaken  the  horizontal  buoyancy  stratification. 
That  is,  the  thermohaline  diapycnal  fluxes  remove  the 
energy  stored  in  horizontal  buoyancy  gradients.  The  sit- 
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uation  is  reversed  below  the  ML  base,  where  the  avail- 
able potential  energy  is  believed  to  be  removed  by  adi- 
abatic processes.  The  difference  in  the  two  cases  is  that 
there  are  strong  diabatic  motions  in  the  ML  at  small 
scales,  while  flows  are  mostly  along  isopycnals  in  the 
ocean  interior.  Expressions  like  the  one  in  (7)  might  be 
a first  step  toward  a better  parameterizations  of  diapy- 
cnal  fluxes  in  the  ML. 

We  implemented  the  nonlinear  diffusive  parameteri- 
zation of  heat  and  salt  in  a simple,  idealized  model  of 
the  ML.  Numerical  simulations  showed  that  buoyancy- 
driven  diffusion  is  complemented  and  locally  enhanced 
by  the  mesoscale  stirring  field.  The  strain  in  the  ve- 
locity field  continuously  creates  thermohaline  gradients 
at  small  scales.  Nonlinear  diffusion  selectively  removes 
buoyancy  gradients.  As  a result  at  any  single  time  typi- 
cal temperature  and  salinity  gradients  are  compensated, 
in  agreement  with  the  observations  of  Ferrari  and  Rud- 
nick  (2000).  Our  model  does  not  include  any  feedback 
between  the  velocity  field  and  buoyancy,  as  due  to  ro- 
tation. An  obvious  direction  for  future  research  is  to 
investigate  whether  this  feedback  plays  an  important 
role  in  the  establishment  of  the  temperature-salinity  re- 
lationship of  the  ML. 

Our  discussion  emphasized  the  role  of  diapycnal 
fluxes  in  removing  horizontal  buoyancy  gradients.  A 
question  arises  as  to  what  maintains  the  long-lived 
buoyancy  fronts  observed  at  some  locations  in  the  ML 
(e.g.  Roden , 1975;  Rudnick  and  Luyten , 1996).  First, 
these  fronts  have  horizontal  scales  of  a few  tens  of  kilo- 
meters. The  buoyancy-driven  fluxes,  discussed  above, 
are  mostly  active  at  scales  shorter  than,  say,  10  kilome- 
ters and  are  weaker  at  larger  scales.  Second,  ML  fronts 
are  believed  to  be  maintained  by  surface  forcing  or  by 
convergences  in  the  velocity  field.  When  this  is  the  case, 
the  diapycnal  thermohaline  fluxes  do  not  remove  com- 
pletely the  buoyancy  gradient.  Instead  an  equilibrium 
is  reached  between  nonlinear  diffusion  and  forcing.  The 
result  is  a buoyancy  front  in  which  temperature  and 
salinity  partly  oppose  each  other,  but  do  not  compen- 
sate. Observations  indeed  show  partial  cancellation  of 
the  thermohaline  components  at  many  ML  fronts. 

We  also  derived  a linear  model  that  reproduces  the 
preferential  diffusion  of  horizontal  buoyancy  gradients. 
The  linear  parameterization  produced  realistic  distribu- 
tions of  the  density  ratio  in  a simple,  idealized  model 
of  the  ML.  It  remains  an  open  question  whether  the 
new  parameterization  would  produce  realistic  diapyc- 
nal fluxes  and  water  mass  conversions,  if  implemented 
in  high-resolution  ocean  models  with  bulk  ML. 
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